Poincare oscilator

$$ \begin{equation} \dot{x} = \gamma * (A - r) * x - \omega * y \end{equation} $$ $$ \begin{equation} \dot{y} = \gamma * (A - r) * y + \omega * x \end{equation} $$ $$ \begin{equation} r = \sqrt{x^2 + y^2} \end{equation} $$

import sympy as sp

# Define the variables and parameters
x, y, gamma, A, v = sp.symbols('x y gamma A v', real=True)
r = sp.sqrt(x**2 + y**2)

# Define r in terms of x and y
r = sp.sqrt(x**2 + y**2)

# Define the system of equations
dxdt = gamma * (A - r) * x - 2 * sp.pi * v * y
dydt = gamma * (A - r) * y + 2 * sp.pi * v * x

# Compute the Jacobian matrix
J = sp.Matrix([[dxdt.diff(x), dxdt.diff(y)],
               [dydt.diff(x), dydt.diff(y)]])
J

$\displaystyle \left[\begin{matrix}- \frac{\gamma x^{2}}{\sqrt{x^{2} + y^{2}}} + \gamma \left(A - \sqrt{x^{2} + y^{2}}\right) & - \frac{\gamma x y}{\sqrt{x^{2} + y^{2}}} - 2 \pi v\\- \frac{\gamma x y}{\sqrt{x^{2} + y^{2}}} + 2 \pi v & - \frac{\gamma y^{2}}{\sqrt{x^{2} + y^{2}}} + \gamma \left(A - \sqrt{x^{2} + y^{2}}\right)\end{matrix}\right]$

Sled in lastne vrednosti

$$ J = \begin{bmatrix} -\frac{\gamma x^2}{r} + \gamma (A - r) & -\gamma \frac{x y}{r} - \omega \\ -\gamma \frac{x y}{r} + \omega & -\frac{\gamma y^2}{r} + \gamma (A - r) \end{bmatrix} $$

$$ Tr(J) = 2 \gamma A - 3 \gamma r = \gamma (2 A - 3 r) $$

Fiksne točke

$$ \begin{bmatrix} \dot{x}\\ \dot{y} \end{bmatrix} = \begin{bmatrix} 0\\ 0 \end{bmatrix} = \begin{bmatrix} \gamma * (A - r) * x - \omega * y\\ \gamma * (A - r) * y - \omega * x \end{bmatrix} $$

$$ \omega * y = \gamma * (A - r) * x $$

$$ -\omega * x = \gamma * (A - r) * y $$

Dobimo dve rešitvi

1. $x = y = 0$

2. $ r = A $

Stabilnost

Točka $x_1 = y_1 = 0$ je stabilna, če je $Tr(J) > 0$ in $Det(J) > 0$ $$ Tr(J)|_{x1,y1} = 2*\gamma*A $$ $$ Det(J)|_{x1,y1} = \gamma^2 A^2 + \omega^2 $$

Vidimo lahko, da velja Tr $Tr(J) > 0$ in $Det(J) > 0$. To pomeni, da je točka $x_1 = y_1 = 0$ nestabilna.

Točka $r = A$. $$ Tr(J)|_{x1,y1} = -\gamma*A $$ $$ Det(J)|_{x1,y1} = \omega^2 $$

Vidimo lahko, da velja Tr $Tr(J) > 0$ in $Det(J) > 0$. To pomeni, da je točka $r = A$ nestabilna.

# Substitute r = A (to evaluate on the limit cycle where amplitude is constant)
J_on_limit_cycle = J.subs(x, 0).subs(y,0)
print("\nJacobian on the limit cycle (x = y = 0):")
J_on_limit_cycle

Jacobian on the limit cycle (x = y = 0):

$\displaystyle \left[\begin{matrix}A \gamma & - 2 \pi v\\2 \pi v & \text{NaN}\end{matrix}\right]$

# Substitute r = A (to evaluate on the limit cycle where amplitude is constant)
J_on_limit_cycle = J.subs(r, A)
print("\nJacobian on the limit cycle (r = A):")
J_on_limit_cycle

Jacobian on the limit cycle (r = A):

$\displaystyle \left[\begin{matrix}- \frac{\gamma x^{2}}{A} & - 2 \pi v - \frac{\gamma x y}{A}\\2 \pi v - \frac{\gamma x y}{A} & - \frac{\gamma y^{2}}{A}\end{matrix}\right]$

eigenvalues = J_on_limit_cycle.eigenvals()
print("\nEigenvalues of the Jacobian on the limit cycle (r = A):")
eigenvalues

Eigenvalues of the Jacobian on the limit cycle (r = A):

{-gamma*(x**2 + y**2)/(2*A) - sqrt((-4*pi*A*v + gamma*x**2 + gamma*y**2)*(4*pi*A*v + gamma*x**2 + gamma*y**2))/(2*A): 1,
 -gamma*(x**2 + y**2)/(2*A) + sqrt((-4*pi*A*v + gamma*x**2 + gamma*y**2)*(4*pi*A*v + gamma*x**2 + gamma*y**2))/(2*A): 1}

trace = J_on_limit_cycle.trace()
determinant = J_on_limit_cycle.det()
print("\nTrace:")
trace

Trace:

$\displaystyle - \frac{\gamma x^{2}}{A} - \frac{\gamma y^{2}}{A}$

#Determinant of J
print("\nDeterminant:")
determinant

Determinant:

$\displaystyle 4 \pi^{2} v^{2}$

import numpy as np
import matplotlib.pyplot as plt

# Parameters
A = 1.0  # Example value for A
gamma = 1.0  # Example value for gamma

# Define a grid of x and y values
x_vals = np.linspace(-2, 2, 100)
y_vals = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x_vals, y_vals)

# Calculate the trace at each point on the grid
Trace = gamma*(2*A-3*np.sqrt(X**2 + Y**2))

# Plot the trace using imshow
plt.imshow(Trace, extent=(x_vals.min(), x_vals.max(), y_vals.min(), y_vals.max()), origin='lower', cmap='viridis')
plt.colorbar(label='Trace')
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'Trace of $\left(-\frac{A \gamma x}{2} - \frac{A \gamma y}{2}\right)$')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
class Oscillator:
    def __init__(self, index, gamma, A, v, initial_x=0.0, initial_y=0.0):
        self.index = index
        self.gamma = gamma
        self.A = A
        self.v = v
        self.x = initial_x
        self.y = initial_y

    def compute_radius(self):
        return np.sqrt(self.x**2 + self.y**2)

    def update_dynamics(self, coupling_x, coupling_y, epsilon, dt):
        r = self.compute_radius()

        # Compute derivatives based on the equations
        dxdt = self.gamma * (self.A - r) * self.x - 2 * np.pi * self.v * self.y + epsilon * coupling_x
        dydt = self.gamma * (self.A - r) * self.y + 2 * np.pi * self.v * self.x + epsilon * coupling_y

        # Update the state using Euler's method
        self.x += dxdt * dt
        self.y += dydt * dt

    def get_phase(self):
        return np.arctan2(self.y, self.x)

    def get_state(self):
        return self.x, self.y

# Initialize an oscillator at position (i, j) with specific parameters
gamma = 10.0
for init_x in [0.2,0.5,0.8,1.5,2]:
    osc = Oscillator(index=(0, 0), gamma=gamma, A=1.0, v=0.5, initial_x=init_x, initial_y=0.0)
    dt = 0.01  # Time step

    phase = []
    x_state = []
    y_state = []
    t_array = np.arange(0,10,dt)
    for time in t_array :
        osc.update_dynamics(0, 0, 0, dt)

        phase.append(osc.get_phase())
        x_state.append(osc.get_state()[0])
        y_state.append(osc.get_state()[1])

    plt.scatter(x_state[0],y_state[0])
    plt.plot(x_state ,y_state,lw=1)
    
# Calculate the trace at each point on the grid
Trace = gamma*(2*A-3*np.sqrt(X**2 + Y**2))

# Plot the trace using imshow
plt.imshow(Trace, extent=(x_vals.min(), x_vals.max(), y_vals.min(), y_vals.max()), 
           origin='lower', cmap='viridis')
plt.colorbar(label='Trace')
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'Trace of $\left(-\frac{A \gamma x}{2} - \frac{A \gamma y}{2}\right)$')
plt.show()
plt.show()

# Initialize an oscillator at position (i, j) with specific parameters
gamma = 1.0
for init_x in [0.2,0.5,0.8,1.5,2]:
    osc = Oscillator(index=(0, 0), gamma=gamma, A=1.0, v=0.5, initial_x=init_x, initial_y=0.0)
    dt = 0.01  # Time step

    phase = []
    x_state = []
    y_state = []
    t_array = np.arange(0,10,dt)
    for time in t_array :
        osc.update_dynamics(0, 0, 0, dt)

        phase.append(osc.get_phase())
        x_state.append(osc.get_state()[0])
        y_state.append(osc.get_state()[1])

    plt.scatter(x_state[0],y_state[0])
    plt.plot(x_state ,y_state,lw=1)
    
# Calculate the trace at each point on the grid
Trace = gamma*(2*A-3*np.sqrt(X**2 + Y**2))

# Plot the trace using imshow
plt.imshow(Trace, extent=(x_vals.min(), x_vals.max(), y_vals.min(), y_vals.max()), 
           origin='lower', cmap='viridis')
plt.colorbar(label='Trace')
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'Trace of $\left(-\frac{A \gamma x}{2} - \frac{A \gamma y}{2}\right)$')
plt.show()
plt.show()

# Initialize an oscillator at position (i, j) with specific parameters
gamma = 0.5
for init_x in [0.2,0.5,0.8,1.5,2]:
    osc = Oscillator(index=(0, 0), gamma=gamma, A=1.0, v=0.5, initial_x=init_x, initial_y=0.0)
    dt = 0.01  # Time step

    phase = []
    x_state = []
    y_state = []
    t_array = np.arange(0,10,dt)
    for time in t_array :
        osc.update_dynamics(0, 0, 0, dt)

        phase.append(osc.get_phase())
        x_state.append(osc.get_state()[0])
        y_state.append(osc.get_state()[1])

    plt.scatter(x_state[0],y_state[0])
    plt.plot(x_state ,y_state,lw=1)
    
# Calculate the trace at each point on the grid
Trace = gamma*(2*A-3*np.sqrt(X**2 + Y**2))

# Plot the trace using imshow
plt.imshow(Trace, extent=(x_vals.min(), x_vals.max(), y_vals.min(), y_vals.max()), 
           origin='lower', cmap='viridis')
plt.colorbar(label='Trace')
plt.xlabel('x')
plt.ylabel('y')
plt.title(r'Trace of $\left(-\frac{A \gamma x}{2} - \frac{A \gamma y}{2}\right)$')
plt.show()
plt.show()

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from tqdm.auto import tqdm as tqdm

# Parameters
N = 40  # Size of the lattice (NxN)
gamma = 1.0  # Self-regulation rate
A = 1.0  # Desired amplitude of oscillations
omega = 1.0  # Natural frequency of oscillators
epsilon = 0.1  # Coupling strength
dt = 0.01  # Time step for integration
steps = 2000  # Number of simulation steps

# Initialize the lattice of oscillators
lattice = np.empty((N, N), dtype=object)
for i in range(N):
    for j in range(N):
        initial_phase = np.random.rand() * 2 * np.pi
        initial_x = np.cos(initial_phase)
        initial_y = np.sin(initial_phase)
        lattice[i, j] = Oscillator(index=(i, j), gamma=gamma, A=A, v=omega, initial_x=initial_x, initial_y=initial_y)

# Simulation function
def simulate_lattice():
    phases = np.zeros((N, N, steps))
    
    for t in tqdm(range(steps),desc='Simulating lattice'):
        # Update each oscillator
        for i in range(N):
            for j in range(N):
                # Calculate coupling from neighbors
                neighbors = [
                    lattice[(i+1) % N, j], lattice[(i-1) % N, j],  # Vertical neighbors
                    lattice[i, (j+1) % N], lattice[i, (j-1) % N]   # Horizontal neighbors
                ]
                coupling_x = sum(neighbor.x - lattice[i, j].x for neighbor in neighbors)
                coupling_y = sum(neighbor.y - lattice[i, j].y for neighbor in neighbors)
                
                # Update oscillator dynamics with coupling
                lattice[i, j].update_dynamics(coupling_x, coupling_y, epsilon, dt)
                
                # Record the phase
                phases[i, j, t] = lattice[i, j].get_phase()
    
    return phases

# Run the simulation
phases_over_time = simulate_lattice()

# Set up the figure and axis for the animation
fig, ax = plt.subplots()
im = ax.imshow(np.sin(phases_over_time[:, :, 0]), cmap="twilight", origin="lower")
plt.colorbar(im, ax=ax, label="Phase (sin)")

# Animation function
def update(frame):
    im.set_array(np.sin(phases_over_time[:, :, frame]))
    ax.set_title(f"Coupled Poincaré Oscillators - Step {frame}")
    return [im]

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=steps, interval=50, blit=True)

# To display the animation inline (if using Jupyter Notebook), use:
# from IPython.display import HTML
# HTML(ani.to_jshtml())

# To save the animation as a .gif or .mp4 file, uncomment one of these:
# ani.save("poincare_oscillators.gif", writer="imagemagick")
# ani.save("poincare_oscillators.mp4", writer="ffmpeg")

plt.show()

Simulating lattice:   0%|          | 0/2000 [00:00<?, ?it/s]

from IPython.display import HTML
HTML(ani.to_jshtml())

Animation size has reached 20977680 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.

Once Loop Reflect

The Kernel crashed while executing code in the current cell or a previous cell. 

Please review the code in the cell(s) to identify a possible cause of the failure. 

Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. 

View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details.

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import matplotlib.animation as animation
N = 20  # Size of the lattice (NxN)
gamma = 1.0  # Self-regulation rate
A = 1.0  # Desired amplitude of oscillations
omega = 1.0  # Natural frequency of oscillators
epsilon = 0.1  # Coupling strength
dt = 0.01  # Time step for integration
steps = 500  # Number of simulation steps

# Initialize the lattice of oscillators
lattice = np.empty((N, N), dtype=object)
for i in range(N):
    for j in range(N):
        initial_phase = np.random.rand() * 2 * np.pi
        initial_x = np.cos(initial_phase)
        initial_y = np.sin(initial_phase)
        lattice[i, j] = Oscillator(index=(i, j), gamma=gamma, A=A, v=omega, initial_x=initial_x, initial_y=initial_y)

# Simulation function
def simulate_lattice():
    phases = np.zeros((N, N, steps))
    
    for t in range(steps):
        # Update each oscillator
        for i in range(N):
            for j in range(N):
                # Calculate coupling from neighbors
                neighbors = [
                    lattice[(i+1) % N, j], lattice[(i-1) % N, j],  # Vertical neighbors
                    lattice[i, (j+1) % N], lattice[i, (j-1) % N]   # Horizontal neighbors
                ]
                coupling_x = sum(neighbor.x - lattice[i, j].x for neighbor in neighbors)
                coupling_y = sum(neighbor.y - lattice[i, j].y for neighbor in neighbors)
                
                # Update oscillator dynamics with coupling
                lattice[i, j].update_dynamics(coupling_x, coupling_y, epsilon, dt)
                
                # Record the phase
                phases[i, j, t] = lattice[i, j].get_phase()
    
    return phases

# Run the simulation
phases_over_time = simulate_lattice()

# Set up the figure and axis for the animation
fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
scatter = ax.scatter([], [], s=30, c='blue', alpha=0.7)  # Scatter points for phases
arrow = ax.annotate('', xy=(0, 0), xytext=(0, 0), arrowprops=dict(facecolor='red', shrink=0.05, width=2))

# Limit the plot to the unit circle
ax.set_ylim(0, 1)

# Animation update function
def update(frame):
    phases = phases_over_time[:, :, frame].ravel()
    x = np.cos(phases)
    y = np.sin(phases)
    
    # Update scatter points in polar coordinates
    scatter.set_offsets(np.c_[phases, np.ones_like(phases)])  # (angle, radius)

    # Calculate the order parameter (mean phase vector)
    order_param = np.mean(np.exp(1j * phases))  # Mean of complex phases
    order_angle = np.angle(order_param)         # Mean phase angle
    order_magnitude = np.abs(order_param)       # Coherence measure (0 to 1)

    # Update the arrow (order parameter vector)
    arrow.set_position((order_angle, order_magnitude))
    arrow.xy = (order_angle, order_magnitude)
    
    ax.set_title(f"Coupled Poincaré Oscillators - Step {frame}")
    return [scatter, arrow]

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=steps, interval=50, blit=True)

# To display the animation inline (if using Jupyter Notebook), use:
# from IPython.display import HTML
# HTML(ani.to_jshtml())

# To save the animation as a .gif or .mp4 file, uncomment one of these:
# ani.save("poincare_oscillators_polar.gif", writer="imagemagick")
# ani.save("poincare_oscillators_polar.mp4", writer="ffmpeg")

plt.show()

# To display the animation inline (if using Jupyter Notebook), use:
from IPython.display import HTML
HTML(ani.to_jshtml())

Animation size has reached 21024736 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.

Once Loop Reflect

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from ipywidgets import interact, FloatSlider
from IPython.display import HTML

class Oscillator:
    def __init__(self, index, gamma, A, v, initial_x=0.0, initial_y=0.0):
        self.index = index
        self.gamma = gamma
        self.A = A
        self.v = v
        self.x = initial_x
        self.y = initial_y

    def compute_radius(self):
        return np.sqrt(self.x**2 + self.y**2)

    def update_dynamics(self, coupling_x, coupling_y, epsilon, dt):
        r = self.compute_radius()
        dxdt = self.gamma * (self.A - r) * self.x - 2 * np.pi * self.v * self.y + epsilon * coupling_x
        dydt = self.gamma * (self.A - r) * self.y + 2 * np.pi * self.v * self.x + epsilon * coupling_y
        self.x += dxdt * dt
        self.y += dydt * dt

    def get_phase(self):
        return np.arctan2(self.y, self.x)

    def get_state(self):
        return self.x, self.y

# Parameters
N = 20  # Size of the lattice (NxN)
gamma = 1.0  # Initial self-regulation rate
A = 1.0  # Initial desired amplitude of oscillations
omega = 1.0  # Initial natural frequency of oscillators
epsilon = 0.1  # Initial coupling strength
dt = 0.01  # Time step for integration
steps = 500  # Number of simulation steps

# Initialize the lattice of oscillators
lattice = np.empty((N, N), dtype=object)
for i in range(N):
    for j in range(N):
        initial_phase = np.random.rand() * 2 * np.pi
        initial_x = np.cos(initial_phase)
        initial_y = np.sin(initial_phase)
        lattice[i, j] = Oscillator(index=(i, j), gamma=gamma, A=A, v=omega, initial_x=initial_x, initial_y=initial_y)

# Simulation function for a single time step
def simulate_step():
    for i in range(N):
        for j in range(N):
            # Calculate coupling from neighbors
            neighbors = [
                lattice[(i+1) % N, j], lattice[(i-1) % N, j],  # Vertical neighbors
                lattice[i, (j+1) % N], lattice[i, (j-1) % N]   # Horizontal neighbors
            ]
            coupling_x = sum(neighbor.x - lattice[i, j].x for neighbor in neighbors)
            coupling_y = sum(neighbor.y - lattice[i, j].y for neighbor in neighbors)
            
            # Update oscillator dynamics with coupling
            lattice[i, j].update_dynamics(coupling_x, coupling_y, epsilon, dt)
    
    # Return the current phases for visualization
    return np.array([[osc.get_phase() for osc in row] for row in lattice])

# Set up the figure and axis for the animation
fig, ax = plt.subplots()
im = ax.imshow(np.sin(simulate_step()), cmap="twilight", origin="lower")
plt.colorbar(im, ax=ax, label="Phase (sin)")

# Animation update function
def update(frame):
    phases = simulate_step()  # Run one step of the simulation with current parameters
    im.set_array(np.sin(phases))  # Update the plot with the new phase values
    ax.set_title(f"Coupled Poincaré Oscillators - Step {frame}")
    return [im]

# Create the animation
ani = animation.FuncAnimation(fig, update, frames=steps, interval=50, blit=True)

# Display the animation inline in Jupyter
HTML(ani.to_jshtml())

# Interactive sliders using ipywidgets for real-time parameter adjustment
@interact(gamma=FloatSlider(value=gamma, min=0.0, max=2.0, step=0.1),
          A=FloatSlider(value=A, min=0.5, max=2.0, step=0.1),
          omega=FloatSlider(value=omega, min=0.5, max=2.0, step=0.1),
          epsilon=FloatSlider(value=epsilon, min=0.0, max=1.0, step=0.05))
def update_params(gamma, A, omega, epsilon):
    # Update the global variables for gamma, A, omega, and epsilon
    for row in lattice:
        for osc in row:
            osc.gamma = gamma
            osc.A = A
            osc.v = omega
    # This function doesn't return anything but will refresh the parameters.

interactive(children=(FloatSlider(value=1.0, description='gamma', max=2.0), FloatSlider(value=1.0, description…